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ABSTRACT 

Emission from high-J CO lines in galaxies has long been proposed as a tracer of X- 
ray dominated regions (XDRs) produced by AGN. Of particular interest is the question 
of whether the obscuring torus, which is required by AGN unification models, can be 
observed via high-J CO cooling lines. 

Here we report on the analysis of a deep Herschel-PACS observation of an extremely 
high J CO transition (40-39) in the Seyfert 2 galaxy NGC 1068. The line was not de¬ 
tected, with a derived 3a upper limit of 2 x 10“^'^ W m“^. We apply an XDR model in or¬ 
der to investigate whether the upper limit constrains the properties of a molecular torus 
in NGC 1068. The XDR model predicts the CO Spectral Line Energy Distributions for 
various gas densities and illuminating X-ray fluxes. In our model, the 00(40-39) upper 
limit is matched by gas with densities ~ 10® — 10^cm“^, located at 1.6 — 5 pc from 
the AGN, with column densities of at least 10^®cm“^. At such high column densities, 
however, dust absorbs most of the 00(40-39) line emission at A = 65.69 p.m. Therefore, 
even if NGC 1068 has a molecular torus which radiates in the 00(40-39) line, the dust 
can attenuate the line emission to below the PACS detection limit. The upper limit is 
thus consistent with the existence of a molecular torus in NGC 1068. In general, we 
expect that the 00(40-39) is observable in only a few AGN nuclei (if at all), because 
of the required high gas column density, and absorption by dust. 
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1. Introduction 


Low-J rotational lines (up to J ~ 7) of CO have long been used as a tracer for molecular 
gas. With the Herschel Space Telescop^ ( [Pilbratt et al. 2010) it was possible to observe CO lines 
of higher J levels, extending the CO Spectral Line Energy Distribution (SLED) to J = 50 (for 
J —7- J — 1). The PACS instrument (Poglitsch et al. 2010) on board Herschel samples the highest 


frequencies, with 13 < J < 50, which we will call the high-J CO lines. These lines with wavelengths 
186.0 pm < A < 52.85 pm arise from states ~ 500 — 7, 000 K above ground, with critical densities of 
~ 10® — 10® cm“®, and are thus good tracers of the Photon Dominated Regions (PDR) in the warm 
and dense Interstellar Matter (ISM) around stars, in shocks, or in the X-ray Dominated Regions 
(XDRs) of AGN. 

AGN unification models require a structure that obscures the UV and optical emission from 
the Broad Line Region of Type 2 AGN from our view (Antonucci &: Miller|1985). Further, compton 


thick AGN, like NGC 1068 (Elvis &: Lawrence 1988), require an X-ray absorbing medium with a 


column density of at least Nh = 10"^^cm“"^, potentially in the form of molecular gas. These two 
effects are often assumed to be caused by a dusty torus, although the gas and the dust are not 
necessarily distributed in the same way. In this paper, we investigate whether high-J CO lines 
can trace the X-ray absorbing gas (from now on called ’the torus’), be it in the form of a classical 
torus (Krolik & Begelman 1988[ [Krolik Sz Lepp 1989) or of a more clumpy and more extended 
cloud structure. The CO emission from a parsec-scale torus was first modeled by Krolik & Lepp 
( |1989 ) (KL89 hereafter), who found that the line fluxes should scale as J® up to Jw 58 (in this 
particular model for NGC 1068), and with absolute fluxes proportional to the total absorbed hard 
X-ray luminosity. More recent torus models have envisioned a dynamical structure of small clouds 


distributed over several parsecs surrounding the AGN (Nenkova et al. 2002 Honig et al. 2006). 


However, these newer clumpy torus models do not yet incorporate a gas phase and have been used 
so far only to predict the dust emission from an AGN. 


In Hailey-Dunsheath et al. (2012) (hereafter HD12) we carried out a detailed analysis of the 


full CO SLED of NGC 1068 up to J = 30 (for higher J only upper limits were found). We concluded 
that two components are responsible for the CO SLED between 10 < J < 31: the ~200 pc diameter 


ring and the northern streamer, both traced by H 2 1-0 S(l) observations (Muller Sanchez et al. 


^Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator 
consortia and with important participation from NASA. 



































2009). Would the molecular torus, if it exists, be observable in the J > 30 lines? In order to answer 
this question, we carried out an additional, deeper, observation of CO(40-39). 


The paper is organized as follows: a description of the existing multiwavelength observations 
of NGC 1068 is given in section [2j Section covers the PACS observations and data reduction and 
presents the full (55 to 200 |J.m) spectrum of the central spatial pixel (spaxel). Section describes 
updates on the PDR-, XDR-, and shock-fits presented in HD12, including the new upper limit on 
CO(40-39) and a small correction to the CO SLED. A brief description of the XDR code used to 
model the torus emission is given in section A discussion of the results follows in section 


2. NGC 1068 


NGC 1068 is a Seyfert 2 galaxy at a distance of 14.4 Mpc, implying a 72 pc/arc second scale 
( ]Bland-Hawthorn et al. 1997). Because of its proximity and luminosity, NGC 1068 is a favorable 
object for testing hypotheses on Seyfert galaxies. 

Early X-ray observations showed emission from the AGN nucleus, but only in scattered light 
( ]Elvis &: Lawren"^ 1988) because the AGN is completely obscured from our view by a column of 
Nh = 10^^ cm“^ (Bauer et al. 2014). The intrinsic 1 - 100 keV luminosity, estimated from the 
observed scattered X-rays, lies between 10^^ ergs“^ and 10^'^ergs“^. The largest uncertainty in the 


X-ray luminosity is the fraction of X-rays that is scattered into our line of sight. Pier et al. (1994) 


give a comparison of different methods to determine this fraction, from which we derive an intrinsic 

' ergs” 


luminosity Li_iookeV = 10^^'® ergs and a spectral slope iT ~ ^ (or photon index T = 2). 


Eurther out from the obscuring torus, the AGN is surrounded by a ring of molecular gas with a 
radius of ~ V — (72-360 pc), also called the circumnuclear disk (CND). The CND is observed in 

CO and H 2 ( Schinnerer et al.||2000 Muller Sanchez et al.||2009 ), and observations of SiO, HCN and 


CN suggest an X-ray driven chemistry (Lepp & Dalgarno 1996 Usero et al. 2004 Garcia-Burillo 


et al. 2010) 


Inside this shell, the near-infrared integral field spectrograph SINFONI provided 0.075^' reso¬ 


lution observations of H 2 1-0 S(l), which reveal two streamers within 0".4 (30 pc) (Muller Sanchez 


et al. 2009). The southern streamer lies between us and the AGN center. The SINFONI observa¬ 


tions furthermore reveal the presence of the northern streamer (0".4 or 30 pc north of the AGN 
center), which may fuel the central black hole. 


A ring of H 2 O masers, with a radius of ~ 20 mas (1.5 pc), is centered at the AGN (Gallimore 
et al.|[l9^ ). The water masers trace densities uh > 10^ cm“^, and could well locate the inner edge 
of a molecular torus. The dust in the torus, on the other hand, is best observed in the IR continuum. 
The NGC 1068 nucleus has been observed several times with the mid-infrared interferometer MIDI 


in order to determine the dust distribution. Lopez-Gonzaga et al. (2014) present the most recent 


observations, which resolve the central 1-10 pc of the AGN. The observations are best fit with 3 
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gray bodies with a Gaussian intensity distribution: a ~ 1.4 pc x 0.5 pc component of 700 K, a 
~ 3 pc X 2 pc component of ~ 300 K, and a component offset to the west (by 1.7 pc) and north (by 
5.6 pc) with a projected size of ~ 14pc x 3.5 pc and a temperature of ~ 360 K. 


3. PACS Observations and Data Reduction 


The CO ladder of the NGC 1068 nucleus has been observed up to J = 30 with the Photo 
detector Array Camera and Spectrometer (PACS) on board the Herschel Space Observatory, and 
is described in HD12. Here we present a deeper observation taken around A = 65.94 pm, where 
we expect the CO(40-39) line with a rest wavelength A = 65.69 pm at a system velocity of ulsr = 
1125kms“^. The observations were taken as part of our GT Key Project SHINING, the OBSIDs 
are listed in Table Among the high J CO lines, 00(40-39) is particularly suitable because 
the spectra of bright galaxies have little contamination from other lines around A = 65.69 pm. 
Moreover, the line traces gas with densities > 10®cm“^ and temperatures > 500K (Table [^, and 
distinguishes a highly excited component from the gas in the CND and streamers. 


We used HIPE 11 for the data reduction. The previous data were analyzed with HIPE 6 
(HD12), so we re-reduced some of the J < 30 lines with HIPE II. The high-J lines detected 
in NGC 1068 are all consistent with arising in the central spatial pixel which covers the central ~ 
(9".4)^ or (700 pc)^ of NGC 1068. Therefore, the data reduction discussed here applies to the central 
spaxel only. The spectra were reduced with the standard pipeline for a “chopped line scan and short 


range scan with background normalization”, including the point source correction factor (Poglitsch 


et al. 2010), and setting the up-sample parameter to 1 and the over-sample parameter to 2. We 
found that the CO(30-29) flux should be slightly modified, we now measure (2± 1) x 10“^^ Wm“^, 
instead of (4.2 ± 1.9) x 10“^'^Wm“^. There is one other small modification we make to the error 
bars published in HD12: the error in CO(22-21) is larger than 30% of the total flux. This line is 
observed at the edges of two adjacent spectral scans, and therefore the uncertainty in the flnx is 
larger than in other lines. The modified line flux is (dj'j^^) x 10“^'^Wm“^. 


The CO(40-39) line observations have been reduced with the same pipeline and parameter 
setting as described above. Eigure[^ shows the region of the CO(40-39) line in the rest frame with 
a spectral resolution of 0.04 pm. There is no clear emission at A = 65.69 pm, but there could be 
an absorption feature around A = 65.6 pm. Potential absorbers are ^®OH Tl^i 2 — 113 / 2 f- — ^ at 
65.69 pm, the CO(40-39) line itself, H 20 ^ 331 — 220 7/2 — 5/2 at 65.61 pm, and NH 2 at 65.57 pm. 
Erom our PACS observations of other galaxies at both 65 pm and 120 pm the line flux ratio of ^®OH 
to ^®OH is always found to be eqnal to or greater than 10 (Gonzalez-Alfonso et al. 2012, 2014). 
Because there is no ^®OH absorption visible in Eigure[^at 65.28 pm, ^®OH can not be responsible 
for the absorption feature at A = 65.6 pm. Moreover, if identified with ^®OH or CO(40-39), the 
line would be blue-shifted by 400kms“^ with respect to the host galaxy, while the highest velocity 
features in the AGN nucleus are the water masers, with ~ 300kms“^ 


Gallimore et al. 1996 


Gonzalez-Alfonso et al. 2012, 2014 
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Table 1. Observed wavelengths 


OBSID Blue band Red band 
(|j,m ) ((j.m ) 


1342203129 

1342203128 

1342203127 

1342203126 

1342203125 

1342203124 

1342203123 

1342203122 

1342203121 

1342203120 

1342259623 

1342259624 


65.0-70.5 

58.5- 66.0 
52.0-59.5 

94.6- 98.0 
91.0-95.1 

87.3- 91.5 
83.2-87.7 
78.9-83.7 

74.4- 79.4 
70.0-75.0 

65.4- 66.5 
65.4-66.5 


130.2-141.2 

117.4- 132.0 
103.8-119.4 

188.6- 196.0 

181.4- 190.4 
174.0-183.0 
166.0-176.0 
157.0-167.7 
148.0-159.0 

139.6- 150.0 


Note. — OBSIDs and covered wavelength of the full Herschel-PACS scan of NGC 1068. All 
observations for the full spectrum have been performed on observational day (OD) 460. The last 
two OBSIDs are from the deeper observations on 00(40-39). 


Table 2. Details of the 00(40-39) line 


transition 

rest wavelength 

upper limit 

^crit 

^upper/ Kq 


(|um) 

(10-1^ W m-2) 

(cm“^) 

(K) 

00(40-39) 

65.69 

2 

10 ^ 

4 X 10^ 
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Absorption by NH 2 is unlikely because it should have several features of comparable intensity 
as indicated in Figure and these are not observed. Moreover, the other NH 2 lines in NGC 1068 
are faint. See Figure for the full (55 to 200 (J.m) PACS spectrum from the central spaxel of 
NGC 1068. H 2 O'’' absorption cannot be excluded completely, with the 65.61 p-m line the strongest 
of the three lines in Figure but it would be an exception in an otherwise emission dominated 
spectrum. A detailed analysis of the PACS spectra of NGC 1068 and other objects will be given in 
Fischer et al. (in preparation). 

Therefore, we assume that the feature around A = 65.6 p-m is an artifact, and treat it as noise. 
This results in a noise level of 0.12 Jy, and a 2>a upper limit on CO(40-39) of 2 x 10“^^ Wm“^, for 
a line width of 250kms“^ (HD12). Figure 0 shows the upper limit superposed on the spectrum. 
If the feature is masked out during the fit, the Scr upper limit becomes 1 x 10“^^Wm“^, and the 
absorbed flux is 4 x 10“^^ Wm“^. The feature could have absorbed the blue part of the 00(40-39) 
line, but it is hard to quantify by how much. Because there is no red wing visible, the artifact 
would have to be coincidently just equal to that of a possible 00(40-39) emission line in order to 
mask it. We find this scenario unlikely and use the aforementioned upper limit of 2 x 10“^^ W m“^ 
for our analysis. The CO SLED with the upper limit and the modifications to HD12 is given in 
Figure 


4. PDR, XDR, and shock models 


With the slightly modified CO(30-29) value and the new CO(40-39) upper limit we repeated 
our PDR, XDR and shock hts, with the same models as in HD12. For comparison we briefly 
summarize the main results from HD12 here. The observed CO SLED was analyzed with Large 
Velocity Gradient- (LVG), shock-, PDR-, and XDR-models. Only the PACS data were htted; 
the 3 < J < 14 CO lines which trace the colder gas were analysed by Spinoglio et al. (2012). 


The PACS SLED is best explained with two components: a moderately excited (ME) component 
with Tkin ~ 169 K, ~ 10^'®cm“^ and a total mass Mh 2 ~ 10®’'^ Mq, and a highly excited 
(HE) component with Tkin ~ 571K, nj/j ~ lO^'^cm”^ and a total mass Mh 2 ~ 10^'^ M©. The 
ME component can arise from X-ray- or shock-heated gas in the CND, while the HE component 
can arise from either the CND or the northern streamer. In order to determine whether the HE 
component arises from the CND or the northern streamer, we followed H12 and compared the PACS 
CO line shifts with the SINPONI observations of H 2 1-0 S(l). The PACS HE lines are blueshifted 
by 59 ± 20 km s“^, which best matches the ~ 25 km s“^ blueshift of the northern streamer. 


To fit our new results, we use the XDR grid from Meijerink et al. (2007). This grid contains 


the CO line emission for gas with densities nn = 10^ — 10®'^ cm“^ and X-ray fluxes Fx = 1.6 — 
160 ergcm“^ s“^. In HD12, the HE component was best fit with Fx = 160ergcm“^ s“^, i.e. the 
maximum value in the grid. This time we used an extended grid with Fx = 9 — 510ergcm“^ s“^ 
(R. Meijerink, private communication). 
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Table and Figure show the parameters of the best XDR-, PDR- and shock-fits. Only 
detected lines were used during the ht, but the best fits were in all cases consistent with the upper 
limits of the high-J CO lines. The XDR fit of the HE component is slightly different from HD12, 
due primarily to the lower value for CO(30-29) and to the extension of the grid: the best fit now 
has an input X-ray flux Fx = 510ergcm“^ s“^, which is again the high end point of the grid. We 
tested whether an even higher X-ray input flux would better fit the HE component, using a different 
XDR grid (described in Section 5.1), which contains X-ray fluxes up to Fx = 10® ergcm“^ s“^. For 
a thermal X-ray input spectrum like in the grid from Meijerink, the CO SLEDs from the 2 codes 
agree within 20% for J < 30. In this model an input flux of Fx = 10^ ergcm“^ s“^ could still fit 
the HE component, but is no improvement over the fit with Fx = 510ergcm“^s“^. 


Our shock modeling uses a two component fit: the ME C-shock model is from Elower & Pineau 
Des Eorets (2010), and the HE C-shock model is from Kaufman &: Neufeld (1996). The lower value 
of the CO(30-29) line has a small effect on the best fit of a shock model to the HE component. The 
HE component can now be fit with a 30 km/s velocity, instead of 40 km/s. The effect on the PDR 
modeling is small (using the PDR model from Meijerink et al. ( 2007| )). A PDR cannot be ruled 
out, but it is hard to match the J = 30 line without overproducing the lower J lines. A single PDR 
component only marginally fits the entire SLED (see Figure]^. Because the kinematics of the lines 
suggest two components, we discard this option. The PDR fit to the HE component (J > 19) alone 
is slightly improved with the lower CO(30-29) flux, but not as good as an XDR or shock fit. 


The upper panel in Eigure also shows the expected CO emission from a pc-scale molecular 
torus with Nh = 10^^ cm“^ and fcovfxn = 0.08, as predicted by KL89, where fcov = 0.25 is the 
surface covering factor of the torus and Fx = Fxaa, x 10^^ergcm“^ s“^. This is further discussed in 
section 15.21 


We conclude that the revised value for CO(30-29) and the new upper limit on CO(40-39) do 
not significantly change the results from HD12. They found that the ME component can arise from 
X-ray- or shock-heated gas in the CND, while X-ray heated gas in the northern streamer is the 
preferred explanation for the HE component, based on the fits and the kinematics. 


5. XDR model of the NGC 1068 nucleus 

In the previous section we have attributed the HE component of the high-J CO SLED to the 
northern streamer, at about 30 pc from the nucleus. If this is correct then we have not detected 
CO emission from the putative torus that is required by AGN unification models. In this section 
we want to address the following question: what does the upper limit on CO (40-39) tell us about 
the existence or non-existence of a torus? If there is such a torus, is it reconcilable with the upper 
limit, and what can be learned about its physical parameters (like size, density, distance from the 
center)? 

In the following sections we model the molecular torus as a semi-infinite slab or as a system 











Fig. 1.— Continuum fit and 3a upper limit: A Gauss profile with a flux of 2 x 10“^^ W m“^ and 
a FWHM of 250 km/s illustrates the 3a upper limit of the CO(40-39) line, at the system velocity 
of NGC 1068 and a rest wavelength of 65.69 pm. The feature around 65.6 pm is considered an 
artifact. 


Table 3. Heating mechanisms 



ME 

HE 

Eull 


XDR 

uh = 10^'^® cm“^ 
Fx = 9ergcm“^s“^ 
A ^ (108pc)2 

= 10®'® cm“^ 

Fx = 510 ergcm“^ s“^ 
A Ri (20pc)2 


1.23 

PDR 


Hff = 10®'® cm“^ 

Go = 10® 

Lpuv ~ 3 X 10®Lq 

Tiff = 10® cm“^ 

Go = 10® 
Lpuv ~ 10^®L© 

2.16 

C-shock 

no = 2 X 10^ cm“^ 

V = 20kms“^ 

A ^ (150pc)2 

no = 10® cm“® 
u = 30 kms“^ 

A Ri (19pc)^ 


1.88 


Note. — Parameters for the best fits with XDR-, PDR-, and shock-models. The XDR and PDR 


models are from Meijerink et al. (2007), with the XDR grid extended to Fx = 510ergcm ^s 


Two models were employed for the shock-fit: the ME component is modeled with Flower &: Pineau 


Des Forets (2010), and the HE component with Kaufman & Neufeld (1996) (like in HD12). A is 


the illuminated area. The values apply to the 2 component XDR- and shock- fit, and to the full 
PDR fit. 
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rest wavelength in //m 


Fig. 2.— The full PACS spectrum of the central spaxel. The wavelength of the brightest NH 2 lines 
are indicated in the spectrum, to show that they are very faint (if detected at all). At 119 pm, the 
flux ratio ^®OH/^®OH « 10, so if no is observed at 65.28 pm, the line can not cause 

the absorption at 65.6 pm. 61 to 63 pm is a known region with artifacts in the relative spectral 
response function. 
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J upper 


Fig. 3.— Best XDR- PDR- and shock-fits to the CO SLED. Red points indicate Plateau de Bure 
Interferometer data (J from 1 — 3, with beam size ranging from 0”.6 to 2”.5; Krips et al. (2011)) 
and Spire data (J from 4 — 13, with beam sizes ranging from 17” to 42”; Spinoglio et al. (2012)), 


blue points represent PACS data (J > 13, with beam size ranging from 11” to 14”; HD12). Only 
the PACS observations were used for the fits. Upper panel: A two component XDR fit with an 
extended grid of Meijerink et al. (2007), and expected CO emission from an X-ray heated molecular 
torus with = lO^^crn”^ and fcovfxAA, = 0.08 (KL89). Central panel: PDR fit to the full PACS 
SLED (thick line), and a PDR fit to the J > 19 observations (thin line). Lower panel: a two 
component shock-fit. Model parameters of all fits are discussed in the text and in Table 
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of such slabs, with strong irradiation from X-rays. The intrinsic X-ray luminosity Li_iookeV = 

1043.5 

ergs ^ indicates input fluxes up to Fx = 10® erg cm ^s ^ are needed to account for gas as 
close as ~ 0.5 parsec from the AGN. Different from the previous section (and from HD12) we here 
use our own XDR code as described in Bruderer et al. (2012) (including benchmark tests), with 
updates in Bruderer (2013). This code, originally developed for the analysis of proto-planetary 
disks, was extended to the higher X-ray input fluxes we need here. 


5.1. XDR grid 


In the model. X-rays with photon energies between 1 keV and 100 keV deposit energy in the 
gas, which has two effects: they heat the dust and drive ionization reactions. The code iterates 
on radiative, thermal and chemical processes, until all are balanced. From this balanced situation, 
the gas temperature and CO rotational line emission is written out as a function of distance and 
column density into the slab. 


In a semi-inhnite slab model, the observer is located on the same side as the X-ray source. 
X-rays enter the semi-infinite slab from one side and deposit their energy in the gas. CO is emitted 
in all directions, but only the CO emission towards the observer is calculated. CO can become 
optically thick between the point of emission and the observer (we did not assume the gas to have 
a large velocity gradient), in this case the attenuated CO emission is given. NGC 1068 is however 
a type 2 AGN, so the gas is located between the observer and the X-ray source. To check whether 
the location of the observer introduces a large uncertainty in our results, we extracted the optical 
depths as a function of column density from the XDR code, for all CO lines. With these optical 
depths, the unattenuated CO emission was calculated as a function of column density, as well as 
the attenuated CO emission towards an observer on the opposite side of the slab. Indeed, for some 
combinations of density and X-ray flux, the CO becomes very optically thick and the observed CO 
emission differs by more than an order of magnitude from one side of the slab to the other. For the 
parameters that are relevant for this paper however, the difference is at most 30% in CO line flux 
at column densities < 10^®cm“^. Another potential issue is geometrical dilution, this is discussed 
further in Section [5.31 


The expected CO emission is calculated for a grid of semi-inhnite slabs with densities n G [10^, 
10^'® ... 10®] cm“^ and Fx G [10“'^,10“^ ... 10®] ergcm“^s“^. The input SED is a power law with 
the shape oc , which best matches the observed XDR spectrum of NGC 1068 (Pier et al 


1994). We also ran the grid for a softer intrinsic spectrum (spectral index of -2 instead of -1). In 


this case, the X-ray absorption and CO emission takes place within smaller column densities than 
for the hard X-ray spectrum. 


Line absorption by dust particles is not taken into account, because this depends strongly on 
geometry and viewing angle. Since we do not know the geometry of gas structures in the central 
few parsec of NGC 1068, we rather assume no line absorption by dust. This can be done safely up 
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to a column Nh = 8.5 x lO^^cm where the dust becomes optically thick at 66 p.m ( Draine|[2003 ). 
At higher column densities, the results from our model are thus upper limits to the line fluxes (see 


Section 5.4) 


In all grid points, the line emission is calculated up to a column of Nh = 10^° cm“^ into the 
slab. High-J CO emission arises from a certain combination of density n. X-ray flux Fx, and column 
density, best represented by the effective ionization parameter ^e// • ^ef f is calculated from formula 

with N 22 the column density in units of 10 ^^ ■ 


12 in 


Maloney et al. 


(1996); ^e// = 47r 


niv; 


.#>+1 ! 
22 


' cm 


-2 


(p = [a — l)/ 7 , in which a = — 1 is the spectral index, and 7 = 8/3 is a parameter that reflects 
that we are interested in the high energy X-rays (E > 1 keV). The high-J CO lines (around J = 
40) are well produced for an effective ionization parameter ^e// ~ 0.02. ^e// is calculated for 3 grid 
points with n = 10^ cm“^, and Fx = lO'^, 10® or 10® ergcm“^ s“^, and plotted in Figure]^ together 
with the gas temperature and the H- and CO- abundances. The central grid point produces most 
CO(40-39) emission: at higher Fx, the gas remains hot and the CO abundance low, while at lower 
Fx the gas quickly becomes molecular and cools rapidly. 


Although H 2 is the main collision partner of CO in a molecular gas, atomic hydrogen also 
contributes to collisional excitation of CO if the gas is partly atomic or ionized. Because collision 
rates between H and CO are currently still uncertain, they are approximated by the H 2 collision 
rate, scaled with the reduced mass of the colliding particles. The next Section treats the CO SLED 
that results from our model. 


5.2. An X-ray irradiated slab 


A model of a torus with a smooth plane parallel slab that is intensively radiated with X-rays 
was already given by KL89, who concluded that most of the torus cooling happens in (very) high J 
CO lines (with J up to ~ 57). A certain column density is however required to produce observable 
CO fluxes, because the gas has a stratihed structure: the first layer, through which the X-rays enter 
the gas, is fully ionized with temperatures of a few thousand Kelvin. The CO forms only further 
into the slab, and cools the gas rapidly. 


In the following discussion the illuminated slab has a covering factor 0.25 since 20-25% of all 
AGN are compton thick (Malizia et al. 2009; Burlon et al. 2011) and its illuminating X-ray flux 
depends solely on the intrinsic luminosity of NGC 1068 (10^®'® ergs“^, assuming an isotropically 
radiating source) and the distance towards the AGN. The absolute emitted CO line fluxes of these 
scenarios are subject to large uncertainties (Rollig et al. 2007), especially at high energy levels. 
One of the reasons for the uncertainties is that the main collision partner for CO, H 2 , forms on 
dust grains. The H 2 abundance thus depends strongly on the dust abundance, its composition, 
and its temperature, which are all uncertain parameters. Moreover, the dust and gas temperatures 
decouple in the extreme conditions considered here, with T^ust ~ a few 100 K, and T, 


a few 

1000 K. So unlike in the study of protoplanetary disks, where absolute line fluxes can be given 


gas 
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in 


within a factor of a few ( ]Bruderer et al.|2012 ), we adopt an uncertainty of 1 order of magnitude i 
the modeled line fluxes. 

A general discussion on what part of parameter space can be constrained by the upper limit 
is given in Section But first we analyse a simple slab model with parameters based on the water 
masers described in Gallimore et ah (1996). In this model, the gas has a density n = 10'^cm“^ 
and an impinging X-ray flux of 10^ ergcm“^ s“^, because the water masers imply gas densities of 
> 10^cm“^ at 1.5 pc from the nucleus. These are the same parameters as used by KL89, and 
with an intrinsic X-ray flux Fx = 10^^'® ergs“^, they represent a torus with inner radius of 1.6 
pc (our intrinsic luminosity and spectral range differs slightly from KL89). Our model prediction 
of the 00(40-39) flux is a factor ~ 8 lower than presented in KL89 - with fcovfxM = 0.08 to 
account for a covering factor 0.25 - as can be seen in the right panel of Figure The reason 
is that KL89 assume the high-J CO lines to be optically thick up to J ~ 57, and therefore the 
luminosity of the line Lj oc oc J^. Instead, we calculate the optical depth as function of column 
density into the slab, and find that the 00(40-39) line does not become optically thick within 
Nh = 10^®cm“^. For larger column densities our model converges towards the KL89 prediction 
(but needs Nh > 10^®cm“^ to match it). At Nh = lO^^cm”^, as used in KL89, the 00(40-39) 
flux differs by a factor 100 between the 2 models. It clearly matters at what column density the 
high-J CO lines become optically thick. 

The upper limit corresponds to a column density of Njj = 1.5 x lO^^cm”^, which is still similar 


to the expected obscuring gas column density in NGC 1068 (Bauer et al. 2014). So although the 

,-2 


00(40-39) line traces the gas we are interested in (n 10^ cm Nh 10^^ cm located within 
a few pc of the AON), the observation did not reach the required sensitivity to constrain the gas 
distribution. 


5.3. A clumpy cloud system 


In the above scenario the gas is distributed smoothly. There are however indications that the 
dust in a torus is distributed in clouds, and most recent torus models treat a clumpy, rather than 
a smooth, dust distribution (see Honig et al. (2006) for such a model on NGC 1068). 


Assuming that also the gas could be distributed in clouds, we approximate a clumpy torus 
with 10 plane parallel slabs. The slabs have the same density as in our previous example and the 
layer closest to the source receives an X-ray input flux of 10® ergcm“^ s“^. All other layers lie at 
larger distance from the source, and are obscured by the intermediate layers, as shown in Figure 
The total column density along the line-of-sight to the source is 10^® cm“^, and the illuminated 
surface of the slabs increases with distance to the source, to keep a covering factor of 0.25. The 
total amount of gas has therefore increased with respect to the previous example. 


The first layer receives the highest X-ray flux, while layers further away from the center receive 
lower X-ray fluxes due to geometrical dilution and due to X-ray absorption in the intermediate 














14 - 



Log Nh Log Log N/j 


Fig. 4.— The gas temperature, effective ionization parameter, and CO- and H- abundances are 
plotted as function of column density into the slab. The density is 10^ cm“^ in all the figures. 
Left: Fx = 10^ ergcm“^ s“^, the gas has become molecular around Nh = lO^^crn”^. Middle: 
Fx = 10® ergcm“^ s“^, the gas becomes slowly molecular and produces most CO(40-39) emission 
around Nh = 10^®cm“^, with ^e// ~ 0.02 (see Figure]^. Right: Fx = 10® ergcm”^ s“^, even at 
Nh = 10^®cm“^, there is too little CO to produce observable CO(40-39) emission. 



Fig. 5.— Left: Cumnlative line emission from an X-ray irradiated slab with a density lO'^ cm“® 
and Fx = 10® ergcm“^ s“^. Right: the CO SLED from a 8.3 pc^ sized structure, at a distance 
of 14.4 Mpc, and with Nh = 10^®cm“^, is plotted together with the prediction by KL89 with 
fcovfxii = 0.08. The difference between the 2 models is a factor 8 for CO(40-39). 


Log abundance 
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layers. The geometrical dilution is accounted for, by taking the grid point with the X-ray flux that 
corresponds to the geometrically diluted - but unattenuated - X-rays from the central source. In 
case this flux lies between 2 grid points, the CO SLED is interpolated between these grid points. 
Each layer thus has a cumulative line distribution attached to it (as shown in the left panel in 
figure]^, from which the line emission is extracted between two column densities. Eor example, 
the emission of the hrst layer could thus be extracted from Nh = 0—7- 10^^cm“^, and that of the 
second from Nu = 10^^ —)• 2 x 10^^cm“^. This accounts for the X-ray attenuation by intermediate 
layers because the XDR model calculates the X-ray attenuation as function of column density. 

The 10 slabs have column densities of 10^^ cm“^. The distance between the slabs is determined 
by the volume filling factor ‘hy: = 0.5 implies that a 0.03 pc thick layer is followed by a 0.03 

pc thick gap. Eigure shows the CO SLED with = 1, 0.5, and 0.2. The 3 plots look very 

similar, there is no sign of a lower output flux or change in SLED for lower <I>y. The reason is 
that the increasing layer surface compensates the geometrical dilution (both oc r^). Only when <I>y 
becomes very small (<I>y ~ 0.05), the outermost layers are located so far from the X-ray source that 
the CO SLED changes significantly. We therefore conclude that the volume filling factor is of little 
importance in our case. 


It is also worth noting that the layered model with <I>y = 1 looks very similar to the single 
slab model that was presented in Section 5.2 Geometrical dilution is thus not important for the 
grid point with n = 10^ cm“^, Fx = 10® ergcm“^ s“^ and Njj = 10^®cm“^. The same applies to 
the other grid points that are used in this paper. 


5.4. Line absorption by dust 

Figure showns that column densities of ~ 10^® cm“^ are needed to match the CO(40-39) 
upper limit in NGC 1068. Within such high column densities, dust absorbs significant fractions of 
the line fluxes, and especially affects the high J GO lines. 

For simplicity, we assume that the gas and dust follow the same distribution, and that the 


4 



4 



Fig. 6.— Schematic overview of the torus models. The cartoon on the left represents one slab with 
a column density 10^® cm“^, the cartoon on the right a 3-layer model with a volume filling factor of 
0.5. The column through these 3 slabs is also 10^®cm“^, but the layers have been scaled in height 
to keep a constant covering factor. The density of all the slabs is constant, while the X-ray flux 
decreases with distance from the central source to account for geometrical dilution. 
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dust resembles galactic dust. 

TT \-2 


The dust cross sections are taken from Draine (2003) and follow 


iT[cm^^per H-atom] oc at the relevant wavelengths]^ In this case the position of the observer 
with respect to the slab does matter for the final SLED. The right panel in Figure]^ shows the CO 
SLED when the torus is located between the X-ray source and the observer, taking line absorption 
by dust into account. The modeled CO(40-39) emission decreased by a factor 10. Because we do 
not know the real dust distribution in NGC 1068, we do not include line absorption by dust in the 
remainder of this paper, but keep in mind that it can attenuate the CO(40-39) flux by a factor 10. 


6. Results and Discussion 

A smooth molecular torus (with n = 10^cm“^, at 1.6 pc from a Li_iookeV = lO^^'^ergs”^ 
source, and with column density Nh = 1.5 x 10^^cm“^) could have been observable in CO(40-39) 
within the sensitivity of the PACS observation. But the uncertainty in the model at these high J 
lines is an order of magnitude. Moreover the dust could attenuate the CO(40-39) line flux by a 
factor 10, so this scenario is not strongly constrained by the upper limit. 

Can we exclude anything else? Is there any grid point where the CO emission exceeds the upper 
limit by more than an order of magnitude? Figure shows the cumulative CO(40-39) emission 
as function of column density and At Nh > 10^^ cm“^ the uncertainty rises above 30% for 
a gas located between the X-ray source and the observer. Line absorption by dust is not taken 
into account. The white, broken, line illustrates where the observed upper limit is matched for an 
illuminated area of dvrr^, with r = 5 pc, 1.6 pc or 0.5 pc, as given in Table 

In all cases, the maximum emission arises at column densities > 10^^cm“^, and at —2 < 
logC"^) < —1.5. So our initial analysis with n = 10^cm“^ and Fx = 10^ ergcm“^ s“^ already 
produced CO(40-39) very efficiently, and most grid points will therefore not be constrained by 
the upper limit. Sticking to a total column density of lO^^crn”^, there are four grid points which 
produce observable amounts of CO(40-39) emission: those with n = 10® and Fx = 10®, n = 10^ and 
Fx = 10®, n = 10® and Fx = 10^, and n = 10® cm“® and Fx = 10® ergcm“^ s“^. We investigate 
here what sizes the gas structures with these parameters must have to match the CO(40-39) upper 
limit. The slab with n = 10® cm“® and Fx = 10® ergcm“^ s“^ is however excluded from the 
analysis, because the CO lines become optically thick, and the observed CO emission depends 
strongly on the position of the observer with respect to the slab. 

Assuming a toroidal gas distribution, with Nh = 10^®cm“^, we calculate the surface covering 
factor (1 for a 47rr^ shell covered with gas) that is needed to match the CO(40-39) upper limit. 
The results are given in Table and are interpreted as follows: 


^Actual cross sections for Milky Way, R_V= 3.1 are found on: http://www.astro.princeton.edu/~draine/dust/ 
dustmix.html 
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Fig. 7.— 10 plane parallel slabs approximate a torus with = 1; 0.5, or 0.2. The total column 
density through the slabs is 10^^ cm“^. Because the increasing slab size compensates the geometrical 
dilution, the total SLED looks similar in these three cases. 



Fig. 8.— Line absorption by dust in a slab with n = 10^ cm“^ and Fx = 10^ ergcm“^ s“^, 
calculated for a torus located between the X-ray source and the observer. Left panel: the broken 
lines indicate unattenuated line emission, while the solid lines indicate attenuated line emission. 
Right panel: dust attenuates the 00(40-39) emission by a factor 10. 
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• A gas layer with n = 10®cm“^ at 0.5 pc from the AGN needs a covering factor of 4.5 to 
be observable, which is not physical. Therefore, such a compact component could not be 
observed in CO rotational lines at the current sensitivity, even if it has a column high enough 
to completely block the X-rays from our view. 

• The slabs placed at 1.6 and 5 pc do give reasonable covering factors of respectively 0.4 and 
0.06. 


MIDI observations of the hot dust in the NGC 1068 nucleus are resolved in two central com¬ 
ponents with scales of 0.5 — 3 pc ( Lopez-Gonzaga et al.|[2M4 ), and a third component that is offset. 
If the gas is distributed at scales similar to these central components, the scenario with a radius of 
5 pc is not valid. This leaves only one scenario left under which the CO(40-39) could be matched: 


the one that was presented in Section 5.2 


If the CO(40-39) line were detected at the upper limit value, it would trace gas between 1.6 
and 5 pc from the AGN, with densities ~ 10® — 10^cm“®, which covers at least 6% to 40% of a 
shell around the AGN. Column densities of 10^® cm“^ are required. Only a few AGN satisfy this 
condition, and the detection of 00(40-39) in an XDR would point to very high column densities 


(although this line can be observed at lower column densities in shocks; Herczeg et al. (2012)). 
Currently, there are no instruments which can perform an even deeper observation of this line, but 
it may be possible with future observatories such as the Space Infrared Telescope for Cosmology 
and Astrophysics (SPICA). 


7. Summary and Conclusion 

The deep observation on the CO(40-39) line in the NGC 1068 central region resulted in an 
upper limit of 2 x 10“^^ W m“^, more than 10 times lower than the previous upper limit (HD12). 
The upper limit does not change the previous conclusions, as to what mechanism heats the Highly 
Excited Component in the CO SLED (HD12): all models (shocks, PDR, XDR) that fit the highest 
excited component in the observed SLED, are consistent with the new J = 40 upper limit. 

We extended an XDR-code, used for protoplanetary disks, with densities up to 10® cm“^ and 
input fluxes up to Fx = 10® ergcm“^ s“^ to investigate in how far the upper limit constrains the 
gas distribution around the AGN. The spectral index of the modeled X-ray spectrum is -1, and the 
slabs have column densities of 10^® cm“^. In the semi infinite slab model, the X-ray source and the 
observer are located on the same side of the slab. The observed CO(40-39) flux differs by at most 
30% if the observer is placed at the opposite side of the slab, as long as Njj < 10^®cm“^. The 
difference increases however rapidly for higher column densities. We investigated the CO emission 
from a gas with density 10^ cm“®, Fx = 10® ergcm“^ s“^, which peaks around J = 40. The modeled 
CO(40-39) flux matches the upper limit for Nh = 1.5 x 10^®cm“^. 

A line detected at upper limit value would trace gas with densities 10® — 10^ cm at 
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Log(Column) Log{Column) Log(Column) 

Fig. 9.— Cumulative CO(40-39) emission in ergs“^ cm“^ sr“^. The white broken line indicates 
where an object with area dvrr^, with r = 5 pc, 1.6 pc, or 0.5 pc respectively (see Table|^, matches 
the upper limit of 2 x 10“^'^Wm“^. Taking Nh = lO^^cm”^ and log('^) = —2, the surface 
covering factor that is needed to match the upper limit is given in Table 


Table 4. Matching the CO(40-39) upper limit 


Fx 

(ergcm“^ s“^ ) 

n 

(cm“®) 

Inner radius 

(pc) 

Outer radius 

(pc) 

Covering factor 
(0-1) 

Required emission 
ergg-i cm“^ sr“^ 

10® 

10® 

0.5 

0.53 

4.5 

1.18 

10® 

10^ 

1.6 

1.9 

0.4 

0.118 

10^ 

10® 

5 

8.2 

0.06 

0.0118 


Note. — Surface covering factors required to fit the J = 40 upper limit face value (line absorption 
by dust is not taken into account). The third column in the table lists the inner radius of the torus, 
assuming that Ti_iookeV = 10^^'® ergs“^. The outer radius follows from a column of 10^® cm“^. A 
10® cm“® density slab located at 0.5 pc from the AGN, needs a covering factor > 1 to be observable, 
which is not physical. The scenarios with 10^ cm“® and 10® cm“® could be observable in CO(40-39). 
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1.6 — 5 pc from the AGN, with column densities of at least 10^^ cm“^, and surface covering factors 
of at least 6% — 40%. Especially because of the high column density required to be observable 
within the PACS sensitivity, we expect to see only few AGN with 00(40-39) emission. 

Because tori are thought to be clumpy structures, rather than a smooth medium, we approx¬ 
imated a clumpy distribution with 10 slabs and volume filling factors 1, 0.5, and 0.2. This did 
not change the CO SLED at 00(40-39), because the increasing layer surface compensates the 
geometrical dilution. A difference is only seen when <l>y = 0.05. 

On the other hand, line absorption by dust is a very important effect. Dust attenuates the 
high J CO lines, which makes a difference of a factor 10 in the modeled 00(40-39) line flux at a 
column density Nh = 10^®cm“^. 

In HD 12 we speculated that there might be no molecular torus, based on the non-detection 
of several high-J CO lines. But we did not have an XDR code that addresses all the caveats and 
uncertainties that are described above. With our new model, adapted to the high densities and 
X-ray fluxes in an AGN torus, considering the position of the observer and line absorption by dust, 
we conclude that the non-detection of 00(40-39) is still consistent with the existence of a molecular 
torus. 
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